Band and Mott Insulators and Superconductivity in Honeycomb-Lattice lonic-Hubbard Model 
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Motivated by superconductivity (SC) in layered nitrides, we study an ionic-Hubbard model on a honeycomb 
lattice, which consists of two sublattices with an energy-level offset, by using an optimization variational Monte 
Carlo method. Parameter values for the Coulomb interaction and the energy-level offset, to realize a band 
insulator in non-doped state, are evaluated, and the SC state is investigated in a carrier-doped band insulator. 
It is found that, in a weakly doped band insulator, spin fluctuation remains and the spin-singlet d-wave SC 
state is realized. Present results support an unconventional spin-singlet SC state suggested by the experimental 
observations in layered nitrides. 



I. INTRODUCTION 

Most of the superconductivity with high-Tc occur in lay- 
ered materials, such as cuprates^i low-dimensional organic 
conductors, and Fe-based pnictideSi^ Because the high-Tc 
superconductivity is often induced close to the magnetic or- 
dered phase, spin fluctuation is considered to be closely re- 
lated to the Cooper pairing. Nitride /3-7\/NCl (M =Hf, Zr) 
doped with carriers is recognized as one of the layered su- 
perconductor with relatively high Tc, which is around 15K 
for Zr and 26K for Hfi^ In a series of materials, the elec- 
tron conduction occurs in MN bilayers with honeycomb lat- 
tice structure, and the superconductivity is induced by dop- 
ing of electrons into the bilayers, namely by the alkali-metal 
intercalation!^— In contrast to most of the high-Tc layered 
superconductivities, this nitride is a band insulator without a 
magnetic order. 

Experimental researches on the layered nitride supercon- 
ductors, LiiJ\f NCI, have revealed many interesting features in 
the superconducting (SC) state. Nuclear magnetic resonance 
(NMR)ii and muon spin relaxation measurement a ' ^ have 
suggested the quasi two-dimensional superconductors. Re- 
gardless of the relatively high Tc, the T-linear specific heat 
has been found to be very small (~' ImJ/mol^)!^^ The iso- 
tope effect is wea k'^-'^ and the SC gap ratio 2A/ k-gTc is 
around 5 for Lio.i2ZrNClii. and Lio.48(THF)3HfNCl^ which 
is unexpectedly high in comparison with the BCS theory. 
The NMRi^ and specific heat measurementsiS suggest that 
the SC gap is of an unconventional spin singlet pairing and 
anisotropic. These results show that the SC in Lij;il/NCl is 
not of a conventional BCS type. A theoretical calculation for 
the SC state in Li^jMNCl was done in the two-band honey- 
comb lattice model where the fluctuation exchange (FLEX) 
method was appliedj2i The dH- id'-type SC gap was suggested 
to be stabilized in the spin fluctuation mechanism. 

We study the ionic-Hubbard model on a honeycomb lattice, 
which consists of two sublattices with an energy-level offset, 
with the conducting plane of Lij:7\/NCl in mind. We adopt 
an optimization variational Monte Carlo (VMC) method. Be- 



yond the previous theoretical approaches, this method en- 
ables us to accurately evaluate Mott and band insulators and a 
metal to insulator transition, despite the coupling strength of 
electron-electron interaction. First, in the non-doped system, 
we determine the parameter values for the Coulomb interac- 
tion and the energy-level offset to realize the band insulator. 
Next, in the system corresponding to the doped band insulator, 
we investigate stabilities of the SC states with several pairing 
symmetries and discuss a mechanism of the paring. Spin fluc- 
tuation survives with sufficient magnitude even for large level 
offsets, for which a band insulator is realized at half-filling. 
We find that the SC with the c?-wave pairing symmetry is sta- 
bilized in the doped band insulating system. 

In Sect, ini the ionic-Hubbard model on a honeycomb lat- 
tice and the formulation of the VMC method are introduced. 
In Sect. Uni the ground-state phase diagram in the non-doped 
system is presented, and the parameter region in which the 
band insulator is realized is determined. In Sect. |IV| the nu- 
merical results for the SC state in a carrier-doped band insu- 
lator are shown. Section [V] is devoted to the summary and 
discussion. A part of the present results was briefly reported 
inRef. 21. 



II. MODEL AND METHOD 

As an effective model for Lij^Tl/NCl, we consider a Hub- 
bard model on a single-layer honeycomb lattice consisting 
of alternating gAIh and gNhsites with an energy-level offset. 
The Hamiltonian is given as follows. 
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FIG. 1: (Color online) Two types of honeycomb lattice models; (a) 
Type I and (b) Type II. 
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Here, t is the electron hopping from a site in sublattice A(B) 
to the nearest-neighbor (NN) sites in sublattice B(A), t' is the 
electron hopping from a site in sublattice B to the NN site in 
sublattice B. We define the number operatore iv^ 



nf^ with nf^ = c^^c^^ where A(=A, B) is an index for the 
sublattices. The level offset for the two orbitals is denoted by 
2A. We do not introduce the hopping from a site in sublattice 
A to the NN site in sublattice A. This is supported by the tight- 
binding fit of the energy bands calculated by the first-principle 
calculation^^ 

For the following technical reasons, we adopt the two types 
of the honeycomb lattice model, termed Type 1 and Type II, 
as shown in Fig. [T] Type I is the so-called "brick model", 
in which a staggered order is characterized by the momen- 
tum (tt, tt). This model is suitable to investigate the spin 
and charge correlations near the metal-to-insulator transition. 
Type II is useful to examine the SC pairing symmetries in a 
carrier-doped system because of the trigonal symmetry in the 
lattice structure. Results are independent of the types of the 
model. 

We utilize an optimization VMC method, which can cor- 
rectly treat the local electron correlation in the whole parame- 
ter space spanned by U and A. As a variational wave function, 
we use the Jastrow type, \I> = P<1>, where <f> is an one-body 
(Hartree-Fock) part, and 7-" is a many-body correlation factor 
As for V, in addition to the well-knowm Gutzwiller (onsite) 
factor give by 



'if 



we introduce an intersite correlation factor T'q given by 



^0 = 11(1 - mQ 
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variational parameters, r runs over all the adjacent sites in the 
NN bond directions, and A' = (A, B) for A = (B, A). After 
all, the trial wave function is given by 



(9) 



It is known that, the binding effect of a doubly-occupied site 
(doublon) to an empty site (holon) is indispensable to describe 
the Mott insulator as well as the SC state, appropriatelyi^a^— 
Regarding $, we consider two trial states: (i) a state which 
has an antiferromagnetic (AF) order in the Mott insulator and 
a charge disproportion (CD) between the A and B sublattices 
in the band insulator, denoted by $ins, and (ii) a BCS state 
<i>Bcs- The first trial state, $ins, is given by diagonalizing the 
Hartree-Fock (HF) Hamiltonian: 



where 



and 



^k(T 
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The operator c^^. is the Fourier transform of c^^. The sum- 
mation runs over the reduced two-dimensional Brillouin 
zone. We define 

= -A + AcD + sgn((T)AAF, (13) 



and 



Bk = 2t'Ecos(k-T'), (15) 



where (X^r') represents a summation for the vectors con- 
necting NN sites in the A-B (A-A) sublattices. We intro- 
duce Aaf and Acd as mean fields (MF) for the AF and 
CD states, respectively. The MF electron densities are de- 
fined as {nt) = N + Acd + sgn(cr)AAF and (ng) = 
N — Acd — sgn(cr)AAF, where TV is an averaged electron 
number per site. We note that Aaf controls the staggered spin 
density in an AF ordered state, and Acd controls the charge 
density in the charge-rich A and charge-poor B sites. 

In order to construct the BCS state $bcs> we introduce the 
following HF Hamiltonian: 



Hbcs = Hk + H 



(m) 
pair ' 



where 



= E </'kcr'^'ko-0k<T, 
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TABLE I: Gap functions g(k) for various pairing symmetries. The wave number vector is defined for the honeycomb lattice model of Type II 
(Fig.O. 



Symmetry 


m 


»(k) 


s 


2 


ffBCS 


di 


1 


ffBCs(cOsfca; - COS ky) 


d2 


2 


gBCs[cos(2fca: + ky) — cos{kx + 2ky)] 


d+id 


1 


gees [cos kx + e'^ cos{kx + ky ) + e^'S' cos ky] 


P 


1 


Qbcs sin fc^ 


h 


1 


Qbcs [sin k^ - sm{k^ + ky) + sin ky] 


f2 


2 


gBCs[sin(2fci, + ky) - sin{kx + 2ky) — sin(fc^ — ky)] 
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FIG. 2: (Color online) Schematic gap functions of (a) s, (b) di, (c) d2, (d) d+id, (e) p, (f) /i and (g) /2 waves. The honeycomb lattice model 
of Type II (Fig.[B is adopted. 
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The pairing term 7^^";). is defined by 
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c^ki) + H.c. 
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2g{k)i4^c\^ + cg^c^kj + H.c, (20) 



where m takes 1 for a pair of electrons between the different 
sublattices (A-B), and takes 2 between the same sublattices 
(A-A or B-B), and g(k) is a gap function for the pairing. The 
chemical potential is fixed at a value determined from the 
tight-binding band of Hk, which depends on the optimized 
parameter Acd- In Eq. ( fTSl l. the MF order parameter of AF, 
Aaf, is not considered, because we are mainly interested in 



the SC state arising in the nonmagnetic doped band insulating 
state. 

By diagonalizing Hk, we have the two bands {a = 

h,l) and the eigen operators e^. By representing c^^^' by ej^ 
I 

k' 



and eL, the Hamiltonian is rewritten as 



Hbcs = E ^k^kleL + E 4a^e£^e^ki- (21) 
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|Ak|' + (A-AcD)2 



-^k = £k - Mo, 



where sgn(Q:) = 1 (— 1) for a ~ h and 



=sgn(a)g(k), (22) 

(23) 
(24) 
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e£ =sgn(a)Vl^k|' + (A- Acd)2 
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From this Hamiltonian, the BCS state is given by 

^Bcs-Ef^k^kKU) ' io), 



(26) 



ka 



in which A'c is a fixed electron number, and the BCS coher- 
ence factor is defined as 



(m) 



(27) 
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As for the gap function g(k), we factorize it as 
g{k) ^ 5BCS^k, 



(28) 



where 2;k represents anisotropy of the gap function and .gees 
is a variational parameter to be optimized. We study a vari- 
ety of pairing symmetries; the s, di, (1-2, d+id waves as can- 
didates for the singlet pair, and the p, /i, [2 waves for the 
triplet pair Formulation and depiction of these gap functions 
are summarized in Tableland Fig. [2] respectively. The di-, 
d+id-, p- and /i-wave functions correspond to the A-B pairing 
(to — 1), while the s-, ^2- and /2-wave functions correspond 
to the A-A or B-B pairing (to = 2). The di wave, correspond- 
ing to the d^2_y2 wave known in the cuprate superconductors, 
has nodes along the ky = ±kx directions. The d+id-wave is 
the anisotropic and fully-opened gap function. 

In the numerical calculations, we adopt 2 x 10^-5 x 10'' sam- 
ples in most of the VMC simulations. The fixed-sampling 
method is used to optimize the variational parametersj2^ The 
error in the energy expectation value is of the order of lO^^t. 
Cluster sizes of 96 < Ns < 216, where Ns is the number of 
site, are used with the periodic -boundary condition and the 
antiperiodic -boundary condition. We take Ng = with 
L = 10-14 for Type I, and Ns = 2L^/3 with L = 12 for 
Type 11. 



III. INSULATING STATE IN NON-DOPED SYSTEM 

We determinate the parameter set of U and A, at which a 
band insulator is realized in the non-doped system. In Fig.|3] 
we show a phase diagram for the most stable state of vj/j^jg 
(= P$ins), at half-filling where N,, = iVs. In the case of t' = 
0, the phase boundary between the Mott and band insulating 
states is located above the boundary in the classical limit (t = 
t' ~ 0), plotted by a dotted line ofU = 2 A. This tendency of 
the phase boundary is also given in the Hartree approximation, 
and reflects the fact that energy reduction from the classical 
limit by quantum fluctuation is larger for the band insulator 
than that for the AF Mott insulator. This is because, in the 
AF insulator, electrons on the sublattice A have little chance 
to transfer to the adjacent sites, owing to the energy 2 A + U . 

When t' is introduced (see the case of t' /t = 0.35 in 
Fig E), the phase boundary shifts to the band-insulator side for 
A/i ^ 1, and shifts to the AF Mott-insulator side for A/t < 1. 
The former is explained by the purturbational picture with re- 
spect to t/U and t' /U; the energy reduction in the AF Mott 
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FIG. 3: (Color online) Phase diagram at half-filling for the most sta- 
ble state of ^ins, in the plane of U /t and A/t. The boundaries be- 
tween the AF Mott and band insulating regions for t' /t = and 0.35 
are plotted by circles and triangles, respectively. The boundaries ob- 
tained by Hartree approximation and at classical limit (t — t' = 0) 
are also depicted by crosses and a dotted line, respectively. The sys- 
tem size of A'^s — 196 is used in the Type-I model. 
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FIG. 4: (Color online) Charge structure factor Ar(q) in the direction 
of q = (0, 0)-(7r, 0) for different values of A at half-filling. The 
parameters are chosen to be (7/t = 7 and t' = 0. Filled symbols 
are for A < Ac, and open symbols for A > Ac, where Ac is the 
boundary between the band and Mott insulating phases. The system 
size of A^s = 196 is used in the Type I model. 



insulator for t'/t ^ 0.35 is 2AtH'UA/{U'^ - 4A2)2 larger 
than that for t'/t = in the case of A ~ U/2. The phase- 
boundary shift for A/t < 1 is explained by a suppression of 
the AF order due to the frustration effect caused by t'. 

Insulating properties are directly confirmed by calculating 
the charge structure factor. 
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where A'^,^ = ^ia- known in the variational theory 
that an insulating gap in the charge degree of freedom is indi- 
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FIG. 5: (Color online) The order parameter for the AF Mott insu- 
lating state Maf and that for the CD state i\/cD for (a) t' /t = 
and (b) 0.35 at half-filling. The parameter value is chosen to be 
U/t = 7. The vertical arrows indicate Ac, the boundary between 
the Mott and band insulating phases. In both (a) and (b), the system 
size of iVs = 100-196 (of L = 10-14) is used in the Type I model. 



cated by the quadratic behavior in the q dependence of A^(q) 
in the Hmit of |q| ^ In Fig. H A^(q) at U/t = 7 is 
shown for several values of A in the direction of q = (0, 0)- 
(tt, 0). For all the values of A, N{q) near the F point seems 
quadratic in q. When A increases, the whole magnitude of 
N{q) is abruptly increased across A/t ^ 1.9. As shown later, 
A/t = 1.9 is a boundary between the band and Mott insulat- 
ing phases. A metallic state is not found within the present 
parameter region of U and A in Fig. [3] except for the Dirac 
point at (A,t/) = (0,0). 

The boundary between the AF Mott and band insulating 
regions is determined by the optimized variational parameters 
Aaf and Acd: we identify the AF Mott insulator for Aaf > 
0, and the band insulator for Aaf = and Acd > 0. In 
order to confirm the realization of Mott and band insulators, 
we check the following two order parameters: 



and 



(30) 
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where S: 



zX 



). We introduce i\/AF as an order 



McD for the CD state where {Nf') > 1 and (A':f ) < 1. 
In Fig. |5] Maf and A/cd are shown as functions of A at 
U/t = 7. In both cases of t' /t = and 0.35, Maf has a 
large value for small A and vanishes at a critical value of A 
(Ac = 1.9^). On the other hand, A/cd is always positive. 
This value gradually increases with A and shows large mag- 
nitude for A > Ac, associated with a small anomaly at Ac. 
It is shown that the magnitudes of A/af and A/cd are almost 
independent of the system size L. These results indicate the 
phase transition between the Mott and band insulators at Ac. 



IV. SUPERCONDUCTIVITY IN CARRIER-DOPED 

SYSTEM 

Now we show the results for the SC state in carrier-doped 
band insulators. Through this section, a value of U is fixed 
at 7t. As for the gap function introduced into the trial BCS 
state ^'bcs (= /''I'bcs). we check all pairing symmetries 
in Table. Ufor all values of U and A studied in Fig. [3] and 
obtain that the BCS states with the di- and d+id-wave gaps 
are stabilized relative to the normal state, 5'bcs(3bcs = 0). 
Therefore, we here present only the results for the di- and 
d+id waves. In Fig. |6|a), a phase diagram for the stable 
state of ^Bcs is presented in the plane of A /t and S, where 
S{= Nc/Ng — 1) is the doping concentration of electrons per 
site. 

First, we briefly touch the results at half-filling. The BCS 
states of the di- and d+id waves are stabilized within the 
wave function of "fBCS- However, we note that the AF in- 
sulating state obtained from 'i'ins is much more stable than all 
the BCS states in the region of A/i < 1.9 for t' /t = and 
A/t < 2.65 for t' /t = 0.35. Then, we focus on the carrier- 
doped cases with t' = 0. The rfi-wave state for <' = is 
stabilized near half-filling, namely for 6 ~ 0.04, as shown in 
Fig Uta). We show the condensation energies for the di- and 
d+id-wave states at 6 = 0.04 in Fig.|7];a). This is defined by 

ational energies per site obtained by the BCS wave function 
^^BCS and that with .gscs = 0, respectively. It is shown that 
Ec for the di-wave state has a finite value in the region of 
A < 2.5i, although the magnitude (~ 10~^t) is small rela- 
tive to the condensation energy obtained in the square-lattice 
Hubbard model near the half-filling (~ 10~^i)i2i These re- 
sults show that in the case of t' = 0, the di-wave state sur- 
vives up to around A = 2.5t, which is larger than Ac = 1.9t 
where the AF Mott insulator to the band insulator transition 
occurs at half-filling. On the other hand, the d+id-wa\e state 
appears only in the region of A < l.Qt, namely it is not stabi- 
lized in the doped band insulator A realization of the d + id 
state at A = and t' = reproduces the previous numerical 
results 1^ 

In order to confirm the SC state in the case of t' ~ 0, we 
calculate the pair correlation function defined as 



(32) 



parameter for spin density wave in the AF Mott insulator and where Aj.(R^) is the creation operator of the NN singlet pair 
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FIG. 6: (Color online) Phase diagram in the 5 ~ A/t plane cal- 
culated in the wave function 'I'bcs for (a) t'/t — and (b) 0.35. 
Filled circles and open triangles represent the parameters where the 
di-wave state and both the di- and d+id- wave states are realized, 
respectively. We note that, as for the parameter spaces indicated by 
triangles, energies for the two BCS states are degenerated. The ver- 
tical arrows indicate Ac, the boundary between the Mott and band 
insulating phases. The parametr is chosen to be C//t = 7, and the 
system size of A'^s = 96 is used in the Type II model. 



along the x-direction, defined as, 



A't At 



(33) 



When P(r) is finite in the limit of |r| — s- oo, one can judge 
that an off-diagonal long-range order is realized. Instead of 
P(|r| — > oo), we average P(r) for r's, which are taken on a 
line segments of (0, L/2)-(L, L/2) and (L/2, 0)-(L/2, L) in 
the Type-11 model with L = 12. We note that P{r) is almost 
constant for |r | > 6 in this system. In Fig.lTJb), we present the 
averaged P(r), termed p^™, as a function of A at 5 = 0.04. 
It is shown that p™'^ for the di-wave SC state is larger than 
that for the c?H-id-wave SC state, and has a large magnitude, 
in particular for A < 2.5t. These results indicate that the SC 
state of the di-wave gap is realized in the carrier-doped band 
insulator. 

In order to examine the origin of the SC state, we show 
in Fig. [8|a) the maximum of the spin structure factor S{q) 
and the charge structure factor N{q) as functions of A for 
6 = 0.04 and t'/t = 0. The spin structure factor is defined as 
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FIG. 7: (Color online) (a) The condensation energies Ec and (b) the 
averaged pair correlation function P'^^'^ calculated in ^bcs for the 
di- and d+id- wave gaps as functions of A/t. The parameters are 
chosen to be 5 = 0.04, t' — 0, and U/t = 7. The vertical arrow 
indicates Ac. The system size of A'^s = 96 is used in the Type II 
model. 



(34) 



and qmax exhibits a wave number vector at which S'(q) or 
iV(q) takes their maxima. As shown in the figure, S'(qmax) 
has a large magnitude in the region of small A, while 
-^(qmax) is almost unchanged by changing A. The maxima 
of S'(q) are located at q^ax = (0, ±7r/6) (see Fig. Ob)). Al- 
though S'(qmax) is rapidly reduced as A approaches Ac, it 
remains a certain magnitude even in the region of A > Ac. It 
is noticed that 5(qinax) shows a similar A dependence with 
pave shown in Fig.lTJb). This fact suggests that the spin fluc- 
tuation is closely related to the di-wave SC correlation. 

Introduction of t' causes a change in the above result for 
the SC state. For t'/t = 0.35, a phase diagram in the plane 
of A/t and 6 is presented in Fig. |6jb). The di- and d+id- 
wave states are stabilized only for 6 < 0.04 and A < 1.0, and 
consequently the SC state does not appear in the region of the 
doped band insulator. To examine the effect of t', in Fig.|9] we 
compare the energy bands for t'/t = and 0.35 in the case 
of A/t = 2. We plot ejj introduced in Eq. ( |25] |, where the 
optimized Acd for U/t — 7 and fiQ for 5 ~ 0.04 calculated 
in this band are adopted. When t'/t = 0, the band is almost 
flat along (27r/3,0)- (7r/3,7r/3), and its small curvature at 
k = (27r/3, 0) produces the large density of state at the Fermi 
energy. On the other hand, the curvature for t'/t = 0.35 is 
larger than that for t'/t = 0. We interpret that the introduc- 
tion of t' decreases the density of state at the Fermi level and 
destabilizes the di-wave SC state. 
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FIG. 8: (Color online) (a) The maxima of the spin and charge struc- 
ture factors, S(q) and A^(q), as functions of A/t. We define qmax 
to be the momenta at which the correlation functions take their max- 
ima. The vertical arrow indicates Ac. (b) A contour plot of 5(q) at 
A/t = 2. The parameters are chosen to be 5 = 0.04, t' — 0, and 
U/t = 7. The system size of A'^g — 96 is used in the Type II model. 



V. SUMMARY AND DISCUSSION 

To elucidate the origin of SC in Lia;A'/NCl, we exam- 
ined stability of the SC state in the ionic-Hubbard model on 
a single-layered honeycomb lattice, using the optimization 
VMC method. In the non-doped system, we investigated the 
A-U parameter space in which the band insulator is realized. 
For the parameter region corresponding to the carrier-doped 
band insulator, we studied the stabilities of the SC state with 
a variety of pairing symmetries. We found that, even in the 
carrier-doped band insulator, the spin fluctuation survives with 
sufficient magnitude. As shown in Fig.[8jb), the spin correla- 
tion function S{q} takes its maximum around q = (0, ±7r/6). 
The di-wave type SC state, corresponding to the d^2_j,2-wave 
in the high-Tc cuprate, is stabilized near the half-filling, in 
comparison with the SC states with other symmetries includ- 
ing the d + id wave. This stability of the rfi-wave SC state 
might be related to the spin fluctuation maximized around 
q = (0, ±7r/6), because these wave vectors satisfy the nest- 
ing condition for one of the two Fermi pockets near the half 
filling. 

Our result supports some experimental observations in 
LiajAfNClfi^ii^ii^ First, the SC phase appears in the electron- 
doped band insulator As shown in Fig. [3] the correlation ef- 
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k 



FIG. 9: (Color online) Energy bands - Mo for t'/t = and 0.35 
along k = (0, 0)-(27r/3, 0)-(7r/3, 7r/3). The parameter values are 
chosen to be A/t = 2 and U/t — 7. The optimized value of Acd 
is used in ejj, and fio is determined from this band structure and 
S = 0.04. 



feet shifts the phase boundary toward the AF Mott insulating 
phase. That is to say, there is a parametr space for the band 
insulator which is stabilized by the quantum fluctuation from 
the Mott insulating phase. We suppose that a series of com- 
pound is located in such a region for the band insulator close 
to the phase boundary. The spin fluctuation is suppressed 
due to correlation effect in the band insulating state, but still 
survives. We propose a possibility that remnant spin fluctua- 
tion in the doped band insualtor induces the superconductiv- 
ity. Second point is about the SC state with unconventional 
spin singlet pair The anisotropic gap was suggested by the 
NMR Knight-shif t'^i'^ and the specific heat measurements J2 
From the present VMC calculations, one possibility for the 
gap symmetry is the rfi-wave shown in Fig.|2]which is stabi- 
lized in the lightly doped band insulator close to the phase 
boundary with the Mott insulator (see Fig. |6]l. This is re- 
lated to the spin fluctuation around q = (0,±7r/6), as men- 
tioned above. Some questions caused from the experimental 
observations in Lia;il/NCl are not resolved within the present 
theoretical analyses. One is the doping concentration depen- 
dence of Tc, which is nearly constant lOK) in a region of 
S > 0.10,i2iS2 although the SC is stabilized only for 5 ^ 0.04 
in our result. To solve this question, further factors to stabi- 
lize the SC might be required. Charge fluctuations due to the 
inter-site Coulomb interactions and the interlayer interactions 
in the conducting plane might be candidates. 
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